Kauffman networks with threshold functions 
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We investigate Threshold Random Boolean Networks with K = 2 inputs per node, which are 
equivalent to Kauffman networks, with only part of the canalyzing functions as update functions. 
According to the simplest consideration these networks should be critical but it turns out that 
they show a rich variety of behaviors, including periodic and chaotic oscillations. The results are 
supported by analytical calculations and computer simulations. 



Random Boolean networks (RBN) were introduced by 
S. Kauffman in 1969 0, [3] to model the dynamics of ge- 
netic and metabolic networks but they are also used 
in a social and economic context 0, and for neural 
networks. Although Boolean models represent a strong 
simplification of the far more complex reality, there exist 
several examples where the modelling of a genetic net- 
work by Boolean variables captures correctly the essen- 
tial dynamics of the system HI- F° r this reason, the 
study of RBNs remains an important step on the way 
towards understanding real networks. 

A random Boolean network is a directed graph with 
randomly chosen links between N binary nodes. We de- 
note the state of a node with <jj = ±1 (we call +1 "on" 
and —1 "off"), and the number of inputs per node with 
K. Each node i is assigned at random an update func- 
tion fi. In this paper, we focus on the case K = 2 and 
on threshold functions 



fi = sign yVjCj,,- = sign (si) 



(1) 



where the sum is taken over the two input nodes for node 



and 



-1 for inhibitory connections and cu = 1 



for excitatory connections. (This version of RBN was 
used for instance in Q.) A connection is excitatory with 
probability p+ and inhibitory with probability 1 — p + , 
leading to the following table (Tab. [I]). 
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TABLE I: The four possible update functions for the model 
used in this paper. The input configuration is given in the 
first column, with \ denoting o"j = 1 and J. denoting Oi — — 1. 
The last row gives the probability for each function, the top 
row gives the name of the function according to the Kauffman 
model. 

In agreement with other authors, we define sign(O) = 1. 
Threshold functions are used not only in the context of 



neural networks, but also in models for genetic networks 
[ll&El- These functions represent four of the 12 can- 
alyzing update functions of Kauffman networks. Cana- 
lyzing functions are those non-frozen functions where at 
least one value of a given input can fix the output of a 
node, irrespective of the value of the second input. All 
nodes are updated in parallel according to the rule 

<Ti(t + l) = fi({<Tj(t)}) = /iK(i),^(t)) • (2) 

Node i depends on the nodes j, namely on node i\ and 
ii- 

The configuration of the system a = {<ti , . . . , <7jv} per- 
forms a trajectory in configuration space. As the state 
space is finite and the dynamics is discrete, some states 
will occur again. If a cycle in state space has a set of 
transient states leading to it, it is called an attractor. 

Kauffman classified the dynamics of RBNs according 
to whether it is chaotic or frozen or critical (as described 
in the review llj). In a chaotic network, a perturbation 
at one node propagates on an average to more than one 
node, leading to long attractors. In a frozen network, 
such a perturbation propagates on an average to less than 
one node, and in the thermodynamic limit N oo only 
a finite number of nodes are not constant after a certain 
transient time. Critical networks are at the boundary 
between these two types of behavior, with a perturba- 
tion of one node propagating on an average to one other 
node. Therefore the difference between two almost iden- 
tical initial states increases like a power law in time. The 
number of nodes that are not frozen on all attractors in- 
creases in a critical network as a power law ~ N 2 / 3 of 
the system size N, as was found numerically in [l2| and 
analytically in 1^, 14 1 . 

It is the aim of this paper to show that this classifica- 
tion breaks down for the simple class of RBNs considered 
here. We find a much richer dynamical behavior with not 
only a frozen and a chaotic phase, but also with two types 
of oscillating phases and several critical points. 

Let us first apply the criticality condition in its sim- 
plest version: For all four update functions, the probabil- 
ity that the output changes if one input spin is flipped, 
is 1/2. Since each node is on an average the input to 
two other nodes, a perturbation at one node propagates 
on an average to one other node, and we should expect 



2 



the model to be critical. This is in agreement with the 
finding that K = 2-RBNs that contain only canalyzing 
update functions (but all of them with the same prob- 
ability) are critical [15| . However, this simple argument 
is based on the assumption that all four possible input 
configurations occur equally often, which may be true at 
the beginning of a simulation run, but may be wrong al- 
ready after one timestep. For this reason, Moreira and 
Amaral [l(| argued that the calculation should be per- 
formed such that the input configurations are weighted 
with their frequencies in the stationary state. 

Let us therefore next apply the rule given by Moreira 
and Amaral and let us determine for what values of p+ 
it predicts that the model is frozen, critical, or chaotic. 
We will then see later that the result is still not correct 
for all parameter values. 

We denote with b t the proportion of nodes in state <Ji — 
+1 at time t. In the thermodynamic limit, it changes 
deterministically according to 

ftt+i = 1- [b 2 t (1 - P+f + (1 - hf p 2 + + 

26t(l-&t)p+(l -*>+)]. (3) 

The expression in the square brackets is the probability 
that an input combination leads to Si = —2. In the 
stationary state, we have b t +\ = b t = b with 



4p* + - 2p+ - 1± -^5 - 12p+ 
2(1 - 2p+f 



(4) 



The sign in the numerator has to be chosen such that 
b € [0, 1], therefore only the positive branch remains, see 
Fig. [T] For p + = 1/2, the denominator vanishes, and the 
stationary solution of Eq. is b = 3/4. For p + = 1, 
we have 6 = and 6=1, with the first solution being 
obviously unstable as it is destroyed by one node in the 
state (Ji = +1. The second solution is a stable fixed point 
of the dynamics. For p + — 0, we have b — (—1 + v / 5)/2. 

The mean number of nodes to which a perturbation at 
one node propagates is in the stationary state given by 
27Ti , with 7Ti being the probability that a node changes its 
state when one input is flipped. We obtain it by adding 
up the probabilities for those input configurations which 
allow a transition between an output +1 and —1 and vice 
versa. This is true for half of the input configurations 
leading to Si — (the first 4 terms in the following equa- 
tion) and for all input configurations for which Sj = — 2 
(the last 4 terms): 

TTl = (l-p + )(l-b)(l-p + )b+p + b P+ (l-b) + 
P+ b(l - p+)b + (1 - P+ )(l - b) P+ (l - b) + 

(1 -p+)b(l -p+)b+p+b(l -p+)(l -b) + 
(1 - p+)(l - b) P+ b +p+{\ - b)p+(l - 6) 
=> 7Ti = b + p + - 2bp + (5) 

For p + = 1/2, we obtain tti = 1/2, for p + = 1, we obtain 
7Ti = 0. For p + — 0, we obtain m — b ~ 0.618. We 



therefore conclude that the model is in the frozen phase 
for p + > 1/2, that it is critical for p + — 1/2, and chaotic 
for p + < 1/2. 

The same result is obtained by calculating the station- 
ary value of the Hamming distance between two identi- 
cal network realizations. The Hamming distance is the 
fraction of nodes for two configurations a, a that have 
different values: D = {AN)- 1 Y^Lifa ~ ^f- If we de- 
note with 7r 2 the probability that a node changes its state 
when both inputs are flipped, the time evolution of D is 
given by 



D 



t+i 



2£> t (l- A)7Tl + D t 2 7T2 . 



(6) 



7T2 in the stationary state is obtained by summing all 8 
combinations leading to Sj € {±2}. It can be written as 

7T2 = 1 - 26(l-2p+) 2 + 26 2 (l-2p+) 2 - 2p+ + 2p\ (7) 

and is 1/2 for p + = 1/2 and 1 for p + = 1. If D t is very 
small, we have 
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2A7T1 , 



(8) 



which allows for the growth of a small perturbation if 
2tti > 1 or p + < 1/2, in agreement with our result above. 
The transition from a stationary value D = to a sta- 
tionary value D > occurs at the same point. 
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FIG. 1: The functions b(p+) 
ary value D(p+) vs. p+. 



9+),7T2(p+) and the station- 



Fig. [5] shows D t for different values of p+ and for given 
initial conditions. One can see that D t approaches 
for large times if p+ > 0.5. Furthermore, one can see 
that D t oscillates with period 2 for the smaller values of 
p + . This oscillation is an indication that the dynamics in 
the "chaotic" phase has some structure, which shall be 
investigated in the following. 

Let us therefore have a closer look at the supposedly 
chaotic phase p + < 1/2. We will see that the dynamics 
is not chaotic at all for sufficiently small p + . The con- 
siderations that have lead to our simple phase diagram 
are flawed. The reason is that we have assumed that 6 
becomes for large times stationary for all p+. In order 
to see that this need not be the case, let us first look at 



3 



0.6 
D(t) 



0.4 



p.=0.05 
p + =0.07 
p.=0.09 
p + =0.11 
p.=0.40 
p.=0.50 
p + =0.60 




if 



FIG. 2: The time evolution of the Hamming distance D for 
different values of p+ when b is not stationary. The curves are 
calculated according to Eqs. ((3}, (0 starting from Do = 0.2 
and bo = 0.5. 
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FIG. 3: The map b t vs. bt+i =1-6* for P+ = 0. The fixed 
point b is unstable as depicted by a sample trajectory. 



the situation where p + = 0: We then have b t +i = 1 — b\ . 
This is a one-dimensional map shown in Fig. [3] The fixed 
point is unstable! Instead of having a stationary point 
with a constant proportion of nodes in the two states, 
the system oscillates between a configuration where all 
nodes are switched on and a state where all nodes are 
switched off. This is not chaotic dynamics at all, but 
very stable dynamics. In order to determine the range 
of p + values, for which the fixed point value of b is unsta- 
ble, we performed a linear stability analysis. The ansatz 
bt+i(b + Sbt) = b + Sbt+i leads in linear order in 5b and 
for p + < 1/2 to 



Sbt+i = -2 (b t (1 - 2 P+ f + (1 - 2p + ) P+ ) 



5b, 



= 11-^5- 12p+ + 8 P 2 + j 5b t =■ M ■ 5b t (9) 

In the last step we used Eq. The fixed point is stable 
if the real part of M is smaller than 1, which is the case 



> (3 - V7)/4 = 



Pcb 



0.0886. 



(10) 



Only above this value does the system have a stationary 
state with constant proportions of nodes being "on" and 
"off". 

We finally investigate in more detail the region p+ < 
p c b, where the proportion of "on" and "off" nodes oscil- 
lates with period 2. For p + = 0, every node oscillates 
with period 2, and we have a global attractor of period 
2. This need not necessarily be the case if b oscillates 
with period 2. The attractor could be much larger, while 
the proportion of off and on nodes oscillates still with 
period two. In order to determine for which parameters 
an attractor with period 2 is stable, we performed again 
a linear stability analysis, but now for two time steps to- 
gether. We assume that the system is on an attractor 
of length 2. Let there be every even time step a propor- 
tion x of "on"-nodes and every odd step a proportion y. 
The time evolution is given by Eq. §3§ , but now we com- 
bine two consecutive time steps. We flip one node and 
look how the Hamming distance grows in comparison to 
the undisturbed system after two time steps. The condi- 
tion that information spreading is critical is in the cycle 
with period 2 



7Tl(aO ' 7Tl(y) 
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(11) 



Combining Eq. (|11|) with the time evolution of x and y 
we obtain three equations 

y=l-[x 2 (l- p+f + (1 - xfp\ + 2x(l - x)p+(l -p + )\ 
x=l-[y 2 (l- p+f + (1 - yfp\ + 2y(l - y)p+(l - p+)} 

i = (x+p + - 2xp + )(y+p + - 2yp + ) 

This system can be solved numerically and we obtain a 
critical value p cn = 0.0657. For p + below this value a 
perturbation at one node will die out and all nodes will 
again blink with period 2. Above this value, attractors 
must be longer than 2. 

We checked these analytical predictions by performing 
computer simulations. In order to identify the transi- 
tion at p cn , we measured the median attractor length. 
As shown in Fig. [U we find with increasing network size 
an increasingly sharp transition. Below the transition at 
p cni the proportion of attractors of length 2 converges 
to some nonzero value with increasing system size, in- 
dicating that cycles of length 2 are stable. Above the 
transition, the median attractor length increases more 
and more rapidly with increasing system size, indicating 
a diverging median. Another finding is that attractors 
become again shorter as the critical point p + = 1/2 is 
approached. 

We also evaluated the frequency of phase jumps in b(t) 
on the attractors. Each phase jump is a deviation from 
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median jump frequency vs. p + 



FIG. 4: Numerical verification for the transition p cn and p c 6, 
both marked by vertical arrows. The upper panel shows the 
median attractor length, the lower panel the median jump 
frequency on each attractor candidate in dependence of p+. 
Each data point corresponds to 5000 sample networks of size 
N with fixed p+ and two initial conditions per realization. 
The time evolution is limited to 5000 computational steps for 
both the transient and the attractor length. The median is 
therefore a more sensible observable than the mean value: If a 
attractor candidate cannot be verified because the time limit 
was already reached we can still calculate the median. 



an oscillation with period 2. The result is shown in the 
lower part of Fig. |H With increasing system size, there 
is an increasingly sharp transition at p c b between zero 
phase jumps and a finite proportion of phase jumps. 

We summarize the different types of dynamical behav- 
ior in the following diagram, Fig. 
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FIG. 5: Overview of the dynamic behavior of the model with 
K — 2 in dependence of the parameter p+. 

To conclude, we have shown that the simplest Thresh- 
old Random Boolean Network shows three different types 
of phase transitions and not just the generally expected 
transition between a frozen and a chaotic phase. For pa- 
rameter values p+ < p cni all nodes oscillate stably with 
period two. For p cn < p + < p c j,, the fraction of on-nodes 
oscillate with period two, but attractors are longer. For 
Pcb < P+ < 1/2, the dynamical behavior is chaotic in the 



sense defined by Kauffman. For p + > 1/2, the network 
is in the frozen phase. 

The lesson to be learned from this study is that the 
dynamical behavior of Boolean networks can be much 
richer than expected from simple considerations. There 
is indeed no reason why some model should not also show 
global oscillations with higher periods or period doubling 
cascades in the temporal behavior of b(t). Real genetic 
networks can be expected to have an even richer dynam- 
ical behavior. If the simple classification into "frozen", 
"critical" and "chaotic" networks fails already in the ran- 
dom model presented in this paper, it will be even less 
suitable for real genetic networks, which have attractors 
with very specific properties related to the function of 
the network. A more sophisticated way of describing and 
classifying the dynamical behavior of Boolean networks 
is therefore required. 

This work was supported by the Deutsche Forschungsge- 
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